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I. INTRODUCTION 

In 1995, the Nobel prize in Economics was awarded to Merton and Scholes for the so-called Black-Scholes (BS) 
model first published in 1973 pQ. It was hailed as a milestone in derivative trading, and led to the development of 
elaborate hedging strategies which promise the safe growth of appropriately composed portfolios of financial assets. 
Unfortunately, however, the usefulness of the formula was based on a simplifying assumption that fluctuations of 
assets follow Gaussian distributions. This made it possible to set up a mixture of assets and options that have a 
chance of growing as a safe investment. In fact, this formula was initially quite successful to a number of professional 
speculators. Later, however, it turned out that frivolous trading based only on Black-Scholes model and done without 
a deeper knowledge of market dynamics can lead to dramatic losses and to the formation of speculative bubbles Wr 
0 . The reason for this is that rare events such as drastic falls in the stock markets are much more frequent than 
would be expected from Gaussian distributions. From time to time catastrophic outlier events which should not have 
happened in hundreds of years can occur. Such events have been named ” black-swan” events, after a popular book by 
N. N. Taleb [5]. The existence of such events is a severe obstacle to all hedging attempts. Normally, one considers price 
changes as the result of random steps of a given finite size and demonstrates that these build up to Gaussian random 
walks. This is a consequence of the central limit theorem (CLT). But in general, some of the steps can be extremely 
large, and the combined random walks are of the so-called Levy type. These possess power-like tails which may not 
even have finite variance. They are encountered in many rare events in nature, such as earthquakes, monster waves 
in the ocean, and giant drops in financial markets. Apparently, we need new option price formulas that incorporate 
the possibility of encountering large drops in prices, and are able to compensate for them by a corresponding rise in 
the price of the derivative. 

Many sophisticated models that go beyond Black-Scholes have been introduced in recent decades: among them 
e. g. models based on Levy distributions |6], truncated Levy distributions [7], Multifractal volatility [8], jump pro¬ 
cesses |9 : , and many other approaches. We would like to focus on models that are based on so-called stable distri¬ 
butions. These have found applications in many scientific fields such as multifractal thermodynamics [10;, quantum 
field theory m, evolutionary systems [12], complex dynamical systems [13], etc. These models usually exhibit self¬ 
similarity and power-law behavior in some particular time interval. We also focus on temporal scaling and show 
that models based on double-fractional diffusion (i.e., self-similar scaling in both spatial and temporal variables) can 
successfully simulate situations, when instant fluctuations cause high short-term volatility. In this case we can redis¬ 
tribute the risk for short and long term options options and the volatility of the model can remain unchanged. This 
approach has considerable possibilities for usefulness in further applications in more complex option pricing scenarios. 
By fitting the data of S&P 500 options, we show that for most of the time the optimal value of temporal-scaling is 
very close to classical diffusion based on Levy flight, but for some particular days, e.g., after sudden drops, the risk 
redistribution can be accounted for by time-fractional diffusion process. 
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The paper is divided as follows: in Section |TT] stable distributions are briefly introduced, in Section [TTT] we discuss 
various definitions of fractional derivatives and their relation to stable distributions of the Levy type. Section |IV 


revises previous models of log-Levy option pricing. In Section |V| we introduce double-fractional diffusion equations, 
and discuss their properties. We also show different ways of representing the solutions. Section VI is dedicated to 
the application of double-fractional diffusion to option pricing and the last section is devoted to conclusions and 
perspectives. 


II. STABLE DISTRIBUTIONS 


Stable distributions (also called Levy distributions) constitute an important class of probability distributions. They 
are form-invariant under convolution, implying that the sum of two random variables governed by Levy distributions 
follow such a distribution themselves. Gnedenko and Kolmogorov [14] showed that such distributions are limiting 
distributions of infinite sums of independent random variables with arbitrary distribution, which is the content of 
generalized central limit theorem. Moreover, it is possible to express their probability density function by Fourier 
integrals. In probability theory, the logarithms of the Fourier transforms are referred to characteristic functions. In 
theoretical physics, these are the well-known Hamiltonian operators H{jp). The stable Hamiltonian operators can be 
expressed as a four-parameter set of the following form 


where 


H a ,p;x,o{p) = \n{e ipx ) = ixp - a a \p\ a (1 - i/3sign(p)u>(p, a)) , 


ui(p,a) 


tan(-7ro!/2) if a/1, 
f In \p\ if a = 1. 


( 1 ) 


Each stable hamiltonian corresponds to the stable distribution, which is denoted as L a ,p.^,a(x). The parameters 
acquire following values: a G (0, 2], ft G [—1,1], a > 0, and x G R. Unless specified differently, we consider the 
standard case, i.e., x = 0 and a — 1, and denote the stable Hamiltonian as H a ^(p). 

For a = 2, regardless of /?, the distribution L, 2 ,p{x) is simply the Gaussian distribution. The parameter x has an 
interpretation of a location parameter, and for a > 1 is equal to the expectation value (x). Parameter a is a scaling 
parameter, for a = 2 equal to the standard deviation. Parameter ft is an asymmetry parameter and the parameter a 
is called stability parameter and determines the overall behavior of the probability distribution function. For a < 2 
the Levy distribution L aj p(x) exhibits power-like Levy tails that behave like #-(<*+!) for positive x (see Ref. [15]), so 

L a ,p(x) ~ o:C a ( 1 + ft)x~^ a ^ for x —>> Too (2) 


with a similar behavior for negative x — oo. The only exception is the case of ft = —1, where for a > 1 exhibits 

the distribution the following behavior 


La, — 1 (%) 


1 ( X \ 2 ( q; - 1 ) 1 

2(a-l)\ac a ) 6XP 


-(«-!) 



(<*-!) 


for x +oo 


(3) 


with some constant c a . The proof can be found in [16:. A similar exponential decay holds for the negative tail and 
f3 = 1. For a < 1, the support is bounded to the interval [aj, oo) for /3 = 1, resp. (—oo,x] for (3 = —1. For a < 2, the 
distribution has infinite moments ( \x\ l ) for l > a. 

The two-sided Laplace transform of L a ^ does not exist, unless /3 = 1, so for SR( A) > 0 (see Ref. [16]) the logarithm 
of Laplace image is equal to 

\n{e~ Xx ) = -Xx-X “a" sec ^. (4) 

From the symmetry argument we can deduce that the expectation value ( e TX ) for r > 0 exists only when ft = —1. 
Sometimes, it is advantageous to use an alternative representation of the stable Hamiltonian, which is 

n Mc (p) = He ipx ) = ixp ~ c\p\ a e isi ^ e %, ( 5 ) 

where c and 6 are uniquely determined by the parameters <a, ft and a m- Accessible values of the parameter 0 satisfy 
condition \0\ < min{<a, 2 — a}. The corresponding area in (cq #)-plane is called Feller-Takayasu diamond, and for 
certain values we obtain extremely asymmetric stable distributions (e.g., for a > 1, the value 0 = — <a + 2 corresponds 
to the case ft = 1 and 6 = a — 2 corresponds to the case ft = —1). Further properties of stable distributions, their 
representations and relations between them can be found in Refs. mmm- 
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III. FRACTIONAL CALCULUS 

Fractional calculus generalizes the classical integral and differential calculus to non-integer integration and differ¬ 
entiation. The motivation was originally to find such an operator, for which the second power of the operator would 
be equal to the first derivative. Subsequently, the goal was to generalize the integral and differential calculus for all 
real values. This section presents a basic overview of some possible definitions. 

Let us begin with the fractional integration, which generalizes the well-known Cauchy formula 

nX nXi rX n -1 1 nX 

/ / ■■■ f{x n )dx n ...dx!= - (x- y) n ~ 1 f{y)dy 

Jx 0 Jx 0 Jx 0 !)• Jx 0 

in a natural way from integer values n to non-integer values v by defining 

xJ-xfWj J (x-y) u ~ 1 f(y)dy. (6) 

The fractional integral is a linear operator that obeys the semigroup property 

7"^l n T^2 _ 7"^l+^2 (* 7 \ 

Xo^x u Xo^x — Xo^x • V ' / 

It is easy to show that ^ p= ^ 0 Z^, which is the baseline for the definition of fractional derivative. 


A. Riemann-Liouville derivative 

The relation between fractional integrals and derivatives suggests the introduction the Riemann-Liouville (RL) 
fractional derivative 


*o Kf(x) ■= (x 0 4^[f}) (x) , (8) 

where \v\ denotes the associated ceiling function, i.e., the lowest integer exceeding u. As with ordinary derivatives, 
these derivatives are linear operators, but contrary to the former, they do not satisfy a composition law like (|7|) , i.e., 

° *oK 2 o xo v?. (9) 


In fact, they share only a few of the properties of ordinary derivative operators. On the other hand, some particular 
choices of derivatives, given by special values of xo, recover some of the typical properties of ordinary derivatives. 
One example of such a fractional derivative can be inferred from the derivative of integer power functions x n . For 
arbitrary integers m, n is the m-th derivative of x n equal to 


d m 

dx m 


(n — m)V 


( 10 ) 


By performing a fractional integration, it is easy to show that the fractional derivative of a monomial has the same 
form as in Eq. (10), when xo = 0. This operator will be denoted V v x = qX^. Moreover, this property is valid for 
arbitrary powers and arbitrary derivatives. As an example, we obtain for n = 0 that 


©xl 


r(i - v )' 


(ii) 


Expression © is equal to zero only for v G IN, because of poles in the Gamma function. 

Sometimes, it is advantageous to discuss the formalism at the level of the Laplace images. If we perform the Laplace 
transform of Def. ([8| , we get 


cm(xy,8] 


M 


raw = s a F( S ) [z>r fc "7(*)L=o 


k =0 


( 12 ) 


where represents floor function, i.e., the largest integer not exceeding x. More details can be found in Ref. m- 
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B. Caputo fractional derivative 

One has to notice that the Riemann-Liouville definition of fractional derivatives has some objectionable properties, 
among them non-zero derivative of constant function and unnatural initial conditions. The latter issue for the most 
part limits the application of the above derivative in physical and other real problems, because the fractional initial 
conditions have no physical meaning. These reasons compel us to introduce a different kind of fractional derivative, 
namely the Caputo (C) fractional derivative, which patches up some of the unwanted properties of RL derivatives. It 
is defined as follows 


SKf( x ) :« 


/ 


f H (y) 


r 01 - v) Jx o (x - y) v+1 ~ M 


d y. 


(13) 


The Caputo derivative is more restrictive regarding its domain, because the function / must have at least \v\ 
derivatives. On the other hand, it recovers the desired properties, because *V X 1 = 0, and the Laplace transform 
is equal to 


LC 

cmf(x)-, a ] = mf\(8) = s v F( S ) ( m ) 

k =0 

Now, natural initial conditions are recovered. In case of the Caputo derivative, we can solve the fractional differential 
equation 


* V xf( x ) = X f( x ) 


(15) 


by introduction of the Mittag-Leffler function 


Ev, c( x ) — 


n ^ V{vn + Q 


(16) 


It is easy to see that the function f(x) = E v ^(\x v ) solves Eq. (15). Eventually, Riemann-Liouville and Caputo 
derivatives can be connected via the relation (to be found e.g. in Ref. [20] ) 


\y\ k 
k =0 


(17) 


C. Riesz-Feller fractional derivative 


Equation (15) is a generalization of another well-known property of derivatives, which holds for exponential functions 


exp(Az) = A m exp(Aa^). (18) 


In the case of the Caputo derivative, this formula is generalized in a way that the exponential function is replaced 
by the Mittag-Leffler function. On the other hand, this property can be preserved, when we send the lower bound to 
minus infinity; the operator 


D v J{x):= lim Xa Trj[x) 


(19) 


obeys Eq. (18) even for non-natural derivatives. Such an operator has to be defined on a different space, because the 
fractional operator of this type is convergent for functions that decay faster than \x\~ u for x —» — oo. This operator 
is called the Riesz-Feller derivative and is denoted by D v x . Interestingly, we get the same operator, when we use the 
Caputo derivative approach. In the case of Riesz-Feller derivative it is advantageous to transform into Fourier image 
because it has the following form 


T{D"f{x)-p] = [£>"/] (p) = J dxe ipx J_ 


d y(x - y) " 1 f(y) = {~ip) a f{p). 


( 20 ) 
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The proof can be found also in [21]. Actually, the Riesz-Feller derivative operator is a special case of the so-called 
Levy pseudo-differential operators, for which the fractional Laplacian (— A) u ' 2 is an important example. This class 
of pseudo-differential operators is defined via the Fourier transform through a Hamiltonian 

[ 6D xf](p) = 'H„,o(p)f(p ), ( 21 ) 

where l~Lu,e{p) is the function defined in Eq. <§. Particularly, for extreme values, i.e., for \0\ = minjV, 2 — v} we 
obtain operators proportional to Riesz-Feller fractional derivatives. For the choice 0 = 0, we have the fractional 
Laplacian. These differential operators are closely related to Levy distributions and Levy processes, because they are 
the solutions for the space-fractional diffusion equation 

Jj ~J{x,t) = e D v x f(x,t ). (22) 

IV. OPTION PRICING UNDER LOG-LEVY PROCESS 

In classic Black-Scholes option pricing theory PQ, the price of an option is calculated with the assumption of an 
efficient market. The evolution of an underlying asset is modeled as a log-normal (or geometric Brownian) process 

d S(t) = fjL S S(t)dt + <T S S(t)dW(t ), (23) 

where W(t) is a standard Wiener process, fis is the drift, and as is the volatility of the process. In the following, we 
will be using notation S t = S(t). The assumption of an efficient market implies that there exists another probability 
measure Q, equivalent to the real probability measure P, under which the discounted price is a martingale, so St = 
e~ fls ^ T ~ t \ST\J r t)Q: where T t is a filtration of the stochastic process at time £, and T is the maturity time. For the 
exponential processes, the change of measures is given by Esscher’s transform [22], which states that the Radon- 
Nikodym derivative is 


dIV t 


e x t 

TZx7\T = exp Xt ~ tM$t] 

\e VP 


(24) 


where /xq = ln(exp(Xi)). Particularly, in case of log-normal distribution with the drift equal to interest rate r, the 
process of price evolution has the following form 


S t = S 0 exp 




t + <rW(t) 


(25) 


The price of the option H can therefore be calculated as [9] 


For the European call option, the pricing formula is, with the notation [x] + = max{x,0}, equal to 

i _(rPi 2 

2 to-2 


(26) 


C(S u t) = e~ rT f dx [S t e rT+x - K] 

J R 


J r “ " VZxra* 

By direct differentiation it is easy to show that the option this price fulfills the celebrated Black-Scholes equation 


(27) 


?SM =rC{s , t) - rS ?^-l^ 


d 2 C(S,t) 
dS 2 ' 


(28) 


with the boundary condition C(St,T) ~ [Sr ~ K\ + ■ 

This option pricing model is the one most used worldwide, with many applications, such as estimation of implied 
volatility m ■ On the other hand, a considerable amount of effort has been made during last two decades on the 
development of more advanced evolution models of the underlying assets. As an example, the fractional Brownian 
motion [24] is an elegant model. Furthermore, the complexity of financial markets motivated many authors to introduce 
even more sophisticated option pricing rules, some examples are provided by Refs. [22, 25, 2(y. 

We adopt the approach introduced in Ref. [6] and assume that the price evolution is driven by the log-Levy model 


d S(t) = rS(t)dt + aS(t)dL a ,p(t). 


(29) 
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Contrary to geometric Brownian process, the prices of log-Levy stable process do not possess all moments (Sf). 
What is more, it is not possible to use Esscher transform. The only exception is the case, when /3 = —1, or similarly 
9 = a — 2, because then the two-sided Laplace transform exists, which is expressed in Eq. ©• As a result of the 
exponential decay of the positive tail of the PDF, all moments exist and are finite. Such a model can better describe 
dramatic price drops on the market, which are more frequent than was envisaged by Black-Scholes theory m The 
price process becomes the following time dependence 

S t = S 0 exp [(r + fi) t + crL a _i (£)] , (30) 

where (i = a a sec The corresponding option price, which is given by ( [26] ) , is equal to 

r r p —ipx 

C{St,t) = e~ rT / dx \S t e rT+x - K} + / dA;- e T[ipp+* a H a ,- 1 ( P )] } (31) 

Jr Jr 2?r 

where r = T — t. By introduction of a new variable z = In St and change of integration variable to y = x + rr + z, 
we obtain 

C(z, t) = e~ rT [ d y[e y -K] + [p.jp(*-v) e T\ip(r+n)-n(ip) a ] = f dy [ e w _ K } + g a (z, r\y, 0), (32) 

Jr Jr 27r Jr 

where g a (z,r\y,0) is the Green function (sometimes also called fundamental solution). By further transformations 
£ = z — y + r(r + /i) and g a (£, r) = e rr g a (z , r\y,0) we obtain the equation for g a in the form 

^ ,T) = -n [ a ~ 2 D^g a ] (e, r), (33) 

together with the initial condition g a (£?0) = £(£)• This equation is a fractional Black-Scholes equation for log-prices, 
and for a = 2, we recover the classical diffusion equation. In the Fourier image, the equation has the form of a 
fractional-diffusion equation 


= H a _i(p) 5 a (p,r). (34) 


V. DOUBLE-FRACTIONAL DIFFUSION 

In some recent works, other models that involve fractional time derivatives have been studied [271129] . It has been 
discovered that complex time scaling introduces more general classes of solutions that exhibit interesting phenomena 
and enables one to price options more realistically. The original motivation was provided by fractional Brownian 
motion, originally introduced in [24] . and later applied to option pricing in Refs. [26] 130 ] , The question at stake is 
which particular fractional derivative is the best when generalizing the Black-Scholes equation, driven by a Green 
function g(£,r), obtained as a solution of a double-fractional diffusion equation 

( K dJ + v[ a - 2 D%})g(t;,T)= 0, (35) 

where we have considered two type of the temporal derivatives denoted by the parameter K , namely C d 1 = *V 1 
(Caputo derivative), or RF d 7 = D 7 (Riesz-Feller derivative). The parameter a is the degree of the spatial derivative 
and corresponds to the stability parameter. The parameter 7 is the degree of the temporal fractional derivative 
(corresponding to parameter v in definitions of fractional derivatives). It is called the diffusion speed parameter , 
because, as discussed below, it influences the speed and type of diffusion behavior. In the following sections, we 
compare two classes of double-fractional diffusion equations, particularly equations with Riesz-Feller time derivatives, 
resp. Caputo time derivatives. Both of these equations are examples of a wide class of two-variable pseudo-differential 
equations, which are usually represented via the Laplace-Fourier (LF) image (which means the Fourier image in spatial 
variable and the Laplace image in the temporal variable) as 


a{s)g(p,s) - ao{s)go{p) = b(p)g(p,s ) (36) 

In our case, i.e. when Eq. ( |36| ) represents a LF image of the double-fractional diffusion equation, then b a (p) = iJ Qj _i(p), 
a 7 (s) = s 7 and ao(s) is determined by the type of derivative, for the Caputo derivative clq (s) = s 7-1 , and for Riesz- 
Feller derivative is a RF (s ) = 1. We shall note that the equation requires one initial condition go (p ), which is the 
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Fourier transform of go(0 = #(£, t = 0). In the following, we consider that go(0 = 5(0? which leads to go(p) = 1. 
This is valid for the case, when 7 < 1 (sometimes called slow diffusion). In the following section we focus on that 
case, because there exists a representation which results in a nice interpretation of double fractional diffusion. On the 
other hand, when we consider 1 < 7 < 2 (fast diffusion ), we can also obtain the same form of equation as in Eq. (36), 
we have only to add the second initial condition of the particular form, namely 


dg_ 

dr 


(£,t = 0 ) = 0 . 


(37) 


The question at stake is, whether the solutions of Eq. (36) can be interpreted as probability distributions, i.e. 
whether they are positive. In Ref. m it is shown that for 7 < 1 are solutions positive for all a G ( 0 , 2 ], while in case 
7 > 1 , the parameters have to obey the condition 0 < 7 < a < 2 . 


A. Composition rule for 7 < 1 and smearing kernel 


In the case when 7 < 1, we can derive a special form of the Green function, which gives a nice interpretation of the 
double fractional Green function as a composition of space-fractional Green functions weighted by smearing kernel 
equal to time-fractional Green function. We begin with Eq. (36). The solution to this equation can, in Fourier-Laplace 
image, be expressed as 


g(p, s) = 


ao{s)go{p) 


a(s) - b(p)' 

Under the assumption that Sft(a(s) + b(p)) > 0 , we can apply Schwinger’s formula and obtain 

POO POO 

g(p,s)= dla 0 (s)e~ la -’ {s) g 0 (p)e lboi{p) = dl g s (s, l)g 0 (l,p). 

Jo Jo 

Therefore, the original (£, r)-dependence is given by 


(38) 


(39) 


9(€,t)~ dlg K (r,l)g 0 (l,^). (40) 

Jo 

Because f R g(£,t)d£ = 1 , we obtain that / 0 °° #k(t, l)dl = 1 . One possible interpretation of the variable l is that it 
represents the generic time-parameter of the system and the function gi(r, l) represents the smearing kernel, so that 
the resultant solution is obtained by integration over all solutions for different time-parameters with the weight factor 
given by the smearing kernel. 

After plugging in for b(p) and go(p), we obtain the solution in the form 

pOO POO 

g(p,s)= c \la 0 {s)e~ la{s) e l ' Ha ’~ l{p) = g K (s,l)g a {l,p), (41) 

Jo Jo 

so the solution is given by the superposition of space-fractional-diffusion Green functions for different times 1. The 
smearing kernel gR(s,l) obeys the differential equation 


d_ 

di 


9k(s, l ) 


- a(s)g K (s,l ) 


(42) 


with an initial condition ^(s, 0 ) = ao(s). 

Riesz-Feller kernel: when we assume that the time-derivative operator is equal the Riesz fractional derivative 
D ^ this results in the expected property that for r < 0 is gitfr, l) = 0 [32 . We find that the Laplace image is equal 
to [33]: 


d. 

dl 


9r F ( s ?0 — ~ s7 9r F ( s ^) 


(43) 


with the initial condition g^ F (s, 0) = 1, leading to the solution g FF (s, l) = e _Zs7 , which nothing else than the Laplace 
transform of the fully asymmetric stable distribution with stability parameter 7 , asymmetry parameter equal to 1 
and the support contained in Rq". The function g FF (t,l) is not normalized, so according to [32], we have 


P OO POO P 

/ dlg* F (T,l)= dl 

Jo Jo J E 


p -ipr 

di>~2—< 

R ^ 




■/. 


^ r 7— 1 

L dp 27r ptwV mUt)' 


(44) 









Thus, we end with 


g RF ^ T) =J o d* 0M)^L 74 (^) 5a (£,O 


(45) 


where L 7j i is the 7 -stable Levy asymmetric distribution. 

Caputo kernel: in case of Caputo fractional derivative *V 7 , the Laplace transform of is equal to gx(s : l) 


o7-l p~ls^ 


. According to Ref. [34] , the inverse Laplace transform is equal to 


1 






l 


r'y 


where M u (z) is the M function of Wright type, which is defined by the infinite series 

(-*)" 


M v {z) = Y 


n =0 


n!T (—vn + (1 — v)) ’ 


Interestingly, the connection of M-function to asymmetric Levy-distributions is provided by the relation 

1 / x \ cv 

cM v \cM v ) 


r v+l 


MA — 




which is valid for v G (0,1), c > 0 and x > 0. Hence, we can rewrite the Green function g c (£,t) as 


9 C ^ T )=j o = l 


d l 


J V-h Ll ' x \V-h 




(46) 


(47) 


(48) 


(49) 


We compare the properties of both smearing kernels in Appendix [A] Now let us turn attention to an alternative 
representation of the Green function, which is more computationally tractable and includes 7 > 1 . 

Representation of expectations: The composite representation can be also used in the case of calculating 
expectation values 


<F(X)> df = [ F(09 DF (Lt )dC (50) 

J R 

Because the composition rule can also be formally used for 7 > 1, we can rewrite the previous expression as 

p p pOO poo 

/ F(0g° F ^,T)^ = / F(0 / dlg K (r,l)g 0 (l,0 = / dl g K (r,l) < F(X) > La{l) (51) 

Jr J r J o J o 

where < F(X) >L a (i) is the expectation value under Levy-stable process in pseudo-time l. This is important in the 
calculation of Risk-neutral measure by Esscher’s transform 

r r°° f °o / ol'k \ 

Hdf = ln{exp(Xi)) DF = In J exp(^) J dl g K (l, l)g 0 (l, C = J dl g K (l, l) exp \-l a sec — J . (52) 


B. Mellin-Barnes representation of double-fractional-diffusion Green function 


It is possible to introduce an alternative, computationally effective representation of double-fractional diffusion 
equation which is based on the Mellin-Barnes m transform. This representation is also valid for the case 7 > 1. Let 
us again begin with equation (36). It has the solution equal to 


g(p,s) 


S 7-« 

s 7 — T-L a ,e{p ) ’ 


(53) 


where k depends on the type of the derivative. For Riesz derivative is k = 7 , for Caputo derivative we have k = 1. 
We shall note that because of computational reasons, we assume here only solutions for £ > 0. The solution for 
negative values can be then easily obtained from the relation g a ,< 9 , 7l «(0 = g a ,-e, 7 ,«( — 0- Therefore, we formally leave 
0 undetermined, even if we assume only extremal cases, for which (3 = —1. 
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FIG. 1. Comparison of Green functions for ordinary derivative (7 = 1), Riesz and Caputo derivative for 7 = 0.9 (slow diffusion) 
and 7 = 1.1 (fast diffusion) for a = 2 and a = 1.6 . The Caputo Green function highlights the peak of the distribution, 
while Riesz-Feller Green function has slower decay in tails of the distribution. Note that for 7 > 1, the green function exhibits 
wave-like behavior with two peaks receding in time. 


According to pa cm, the inverse Laplace transform of Eq. ( [53] ) is equal to the Mittag-Leffler function 

ff(p, t) = t k ~ 1 E 1 ^ k {n a ,e(p)r 7 ) . (54) 

The Mittag-Leffler function can be represented through the Mellin integral transform which is defined as 


/* 00 

A d[g(x);s] = / g{x)x s ~ 1 dx 

J 0 


and the inverse transform is defined as (for some c given by the Mellin transform theorem [37] ) 

rC+lOO 


-1 pC- 1-too 

9 (x) = — V/[ 5 ](s)x“ s ds. 

Z7TI J Q — ioQ 


(55) 


(56) 


This representation can provide the Green function in a more tractable form which can be better exploited in numerical 
computations. The Mittag-Leffler function can be expressed as a complex integral 


E, 


1 , c+ . oo r(sW _ s0 s ; 

aAz) - 2ni J c _ ioo T{b — as') ( 


(57) 


where 0 < SR(c) < 1 . Plugging into the equation (54), we find that 


K— 1 /»C+Z OO 


t k ~ 1 r 


W(i - »') 


— fi\p\ a exp i — 


T(k — 7 s') 

Transforming the variable p back to variable £, we obtain 

^-1 r c+ioo r(,s / )r(i - s') 


iirO sign(p)\ 

"2 J 


ds'. 


(58) 


9(£,t) = 


^-1 /*Ct 

■*ICI Jc-i 


2tt*|C| Jc-ioc r(/t - 7 s')r(s'a) [ 


ds'. 


(59) 
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Green functions for a =1.8 




Semi-log plots of GF 





FIG. 2. Green functions (left), semi-log plots of g DF (x) (center) and corresponding option prices (right) for a = 1.8 (top line) 
and a = 1.6 (bottom line) and comparison with the Black-Scholes model (grey dashed lines). For some particular choices 
of parameters, there exist regions, where the option prices are cheaper than the Black-Scholes model and vice versa. The 
parameters are set to a = 1 and r — 1. 


After change of variables as' = s and taking into account the normalization (44), which can be in both cases written 
we finally arrive at the normalized Green function 


as 


9 Ut {i,r) = 


r(«) 


1/ 

iUc 


c+zoo 


r (4) r (i~ f)r(i-s) 

2ani ^ r (« - * 8 ) r (^) r (l - ^) 






d s. 


(60) 


From Eq. (60) it is apparent that #(£, r) = -^j^g ^ 77 ^, 1^, so the Green function has the expected scaling with the 

temporal scaling exponent equal to The ratio D = ^ is called the diffusion scaling exponent. Differences between 
Riesz-Feller and Caputo Green functions for various values of a and 7 are displayed in Fig. [l] 


VI. DOUBLE-FRACTIONAL OPTION PRICING MODEL 


The solution of option-pricing equation driven by double-fractional diffusion under interest rate r and dividend 
yield q can be obtained by integrating over all scenarios [St — iF] + , so it reads 


^7O',7,ft) (*^t 5 -^5 ^ 


L 


~ rT ' d y 

R 


g e r O-g+/d+2/ _ j{ 


n + 


9 {v,r)m 


- e 


L 


d y 


S e T ( r - q +n)+ y _ 


i+ r(«; 


If 

iy Jc 


c+zoo 


r(|)r(i-|)r(i- a ) 

2 a ™yJc-ioo r(«-^)r(^i)r(i-^) 




ds. (61) 


In Fig. [ 2 ] Green functions for different pairs (a, 7 ) and semi-log plots are displayed for better tail-behavior illustration. 
Option prices derived from Green functions are also presented. The corresponding put price can be obtained through 
put-call parity relation, which reads: 


P(a, 7,0 K, r) = C (a , 7>#6) (S t , K, r) - S t e~^ + Ke 


(62) 


A. Model calibration for S&P 500 options traded in November 2008 

We calibrate our model on the data of S&P 500 options that were traded during November 2008. The choice of 
this period is mainly because of the financial crisis, which brought about phenomena that are potentially interesting. 
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All options 

parameter 

Black-Scholes Levy stable Double-fractional 

a 

1.493(0.028) 1.503(0.037) 

7 

1.017(0.019) 

cr 

0.1696(0.027) 0.140(0.021) 0.143(0.030) 

AE 

8240(638) 6994(545) 6931(553) 

Call options 

parameter 

Black-Scholes Levy stable Double-fractional 

a 

1.563(0.041) 1.585(0.038) 

7 

1.034(0.024) 

cr 

0.140(0.021) 0.118(0.026) 0.137(0.020) 

AE 

3882(807) 3610(812) 3550(828) 

Put options 

parameter 

Black-Scholes Levy stable Double-fractional 

a 

1.493(0.031) 1.508(0.036) 

7 

1.047(0.017) 

cr 

0.193(0.039) 0.163(0.034) 0.163(0.037) 

AE 

3741(711) 3114(591) 2968(594) 


TABLE I. Estimated mean values and standard deviations of model parameters ( a , 7 , cr) and mean aggregated error (MAE) 
for three considered models, i.e., Black-Scholes model, Levy stable model with ordinary time derivative and Double-fractional 
model. The statistics are calculated for all options and separately for call options and put options. Apparently, the mean 
value of 7 is very close to one for all options. But when call options and put options are analyzed separately, the 7 -value 
is larger than one. The mean square error of Double-fractional model exhibits significant improvement with respect to the 
Black-Scholes model, but compared with Levy stable model, it shows only a little improvement. This is caused by the fact, 
that the assumption of constant parameters does not fully describe the complex behavior of option markets. The situation 
improves if we calibrate the model for call and put options separately. In Fig. [3] it is shown that on some particular days, the 
improvement of aggregated error is more significant. 


We follow the methodology of Carr and Wu [ 6 ] and try to find such a triplet (a, 7 , a) that minimizes the aggregated 
option price error 


^-^model — ^ ^ l^model (^market) (63) 

rer,Kef c 

of all out-of-the money options, so 

(ao, 70 , cr 0 ) = arg min AE DF (a, 7 , cr). (64) 

(a,7 ,cr) 

We make the optimization for each trading day. We have chosen the out-of-the-money options, because in-the-money 
option prices are more determined by the boundary conditions in option pricing formula, rather than by particular 
underlying diffusion model [ 6 |. The statistics of calibrated parameters is listed in Tab. [I] for the Black-scholes model, 
Levy stable model and Double-fractional model. Because of the fact that 7 was close to 1, the choice of derivative 
type did not have a large impact on the solutions and the results are practically the same, therefore the results 
are presented only for Caputo derivative, i.e., k = 1. The parameter a is close to 1.5 in all cases, only for call 
options it fluctuates around 1.6. The analysis shows that the double-fractional model brings more improvement 
when it is fitted separately for call and put options. This could be connected to the discussion about the general 
validity of put-call parity and market efficiency [38]. Fig. [ 3 ] shows estimated parameters for every day and also the 
ratio between aggregated errors of Levy stable model and Double-fractional model p = ^Edf ' Particularly for put 
options we can observe a noticeable improvement. Additionally, another interesting phenomenon connected with 
the decrease of a parameter can be observed. In such a situation the 7 parameter also decreases, even to values 
smaller than one, which could be interpreted as the risk redistribution from long-term options to short term options. 
The parameter £4, which is the ratio between parameters 7 and <a, expresses the temporal scaling parameter of the 
system. This parameter produces more stable behavior, pointing to the simultaneous changes in both parameters 7 
and a. The differences between option prices in the case of Double-fractional and Black-Scholes model are presented 
in Fig. [4] One can observe that the price differences are relatively small, even in some cases are Double-fractional 
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S lability pa lame ter a Diffusion spee d pa lame tel y 



FIG. 3. Estimated values of stability parameter a , diffusion speed parameter 7 , scaling exponent and the ratio p of aggregated 
errors between Levy model and Double-fractional model for each particular day, done for all options, resp. for calls and puts 
separately. We can observe that in case of calls and puts the benefit of the Double-fractional model is more significant. In the 
case when we observe a drop in a parameter, we can simultaneously observe a drop in values of 7 , which points to the risk 
redistribution from long-term to short-term options. This phenomenon could be an indicator of market regime change. The 
parameter Q measures the ratio between 7 and a and corresponds to the temporal scaling exponent, so g(£,r) ~ r n . We shall 
note that for BS model Q, = |, which corresponds to Brownian motion. The graph shows that Q exhibits more stable behavior, 
especially when estimated on both put and call options simultaneously. 


prices below prices estimated by the Black-Scholes model. This is caused by the fact that the Double-fractional Green 
function redistributes the probability of particular scenarios. In cases, when the option execution is less probable 
than in case of Black-Scholes model, we obtain cheaper option price. Nevertheless, with increasing maturity time the 
Double-fractional prices become more expensive, but the difference is still not dramatic. We discuss the topic of risk 
redistribution and hedging in the next section. 


B. Risk Redistribution and Hedging Policy under Double-fractional Diffusion Model 


We notice a few properties of option prices driven by double-fractional diffusion. Indeed, we recover classical BS 
model for a = 2 and 7 = 1 . When a < 2, the underlying return distribution is skewed and we have polynomial decay 
for the negative tail of the distribution. This results in an adjustment of the option price. A greater probability of 
fall of the underlying asset price results in an increase of option price for K < F [315], and options, for which K > F, 
become cheaper (both puts and calls). Similarly, the parameter 7 plays analogous role in temporal risk redistribution. 
For 7 < 1, options with short expiration period become more expensive, while options with long expiration period 
become slightly cheaper. This behavior can be observed in situations when we face some kind of unexpected or 
sudden change of regime, such as a black day on the market, the bankruptcy of a company trading on the market, 
a natural disaster, etc. Indeed, for options with long expiration long-term equilibrium volatility is more important. 
Nevertheless, for options with short expiration such jumps and short-term uncertainty are the most important factors 
for price estimation. On the other hand, for 7 > 1, the diffusion is faster than in case of space-fractional diffusion, and 
the options with long maturity time. We should note that the change of parameters a or 7 does not change all option 
prices in the same direction; there are always options which become cheaper and more expensive. This essential role 
is played by the parameter <r, which is the volatility of the system. 

The options are usually used for hedging the risk coming from the random nature of price evolution in financial 
markets. In the simplest case, with one underlying asset and one option, we use the <fi-hedging strategy [40 . We 












13 


-- - DF 7=0.25 

200 

-- - DF 7=025 

- - BS 7=0.25 


BS 7=025 

DF 7=0.75 

150 

- - DF 7=0.75 

BS 7=0.25 

BS 7=025 

- DF 7=1.5 


- DF 7=1 5 

- BS 7=1.5 

IOO 

BS 7=1i 


yyj' 


0.12 



0.10 



--- DFt=0.25 --- BSt=0.25 ---- DFt=0.75 


0.10 


BSt=0.75 ---■ DFt=1.5 BSt=1.5 

0.0B 

0.08 







0.06 

0.06 







0.04 

0.04 




0.02 



0.02 



11 

11 

11 

11 

11 

11 

11 

1 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

11 

K 


--- DFK=0.8S 
— BSK=0.8S 
■■■ DFK=S 
■■■ BSK=S 


■ DFK=1.2S 

■ BSK=1.2S 


0.4 0.6 


0.B 1.0 1.2 1.4 


FIG. 4. Top: Estimated call and put prices of S&P 500 calibrated by Double-fractional (DF) model and Black-Scholes (BS) 
model as functions of strike price K for several maturity times r. The market spot price is S — 1000. The model parameters 
are listed in Tab. [I] One can observe that the price differences are relatively small, increasing with increasing r, which is caused 
by different scaling exponents D. In some cases the price of DF option is even cheaper than in case of BS. 

Bottom: Optimal policies 0*(A, £) obtained from estimated parameters as functions of K, resp. r for Double-fractional model 
and Black-Scholes A-hedging. 


create a portfolio of an option 


11(5, t) = C{S , t) - 0(5, t)S(t) (65) 

where 0(£) is the amount of stocks used to hedge the short position. There exist several risk measures. Nevertheless, 
we work with the most popular risk measure defined as variance of portfolio between times to = 0 and T 

n = ((An(5, t)) 2 ) = (([S-K}+- C(S 0 , K, T ) - 0 A S) 2 ) ( 66 ) 

The advantage of using this risk measure is that the optimal policy is easily expressible. The minimal risk is therefore 
determined by 


dR „ 

50 - ' 


In Ref. [41] it was shown that the optimal hedging policy is given by 

1 1 r°° 

= ^((S 0 - S)[S - K}+) = ^ / dS(So-S)[S-K]+g(S,T\S 0 ,t ) 
a a Jo 

where a 2 is volatility of stock 5. For Black-Scholes model we obtain the well-known A-hedging rule 


(67) 


( 68 ) 


<f,% s (S,t) = A(S,t) = ^^-. (69) 

The difference in optimal policies of Black-Scholes model (A-hedging) and Double fractional model is shown in Fig. [4] 
Due to the risk-redistribution, the resulting policies are also changed in order to minimize the risk in the environment 
with double-fractional diffusion. As with the prices, the risk is, in the case of changing parameter a, redistributed in 
the spatial axis (prices) and in the case of parameter 7 in the temporal axis (time to maturity). 
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VII. CONCLUSIONS AND PERSPECTIVES 


A novel method of option pricing was proposed on the basis of fully asymmetric Double-fractional diffusion and a 
special example was calibrated for S&P 500 options traded in November 2008. The presence of power-law behavior in 
prices has been discussed before in case of cotton prices [4Z and options IS 1 E]. There have been various alternative 
attempts to reduce the risks in portfolios. Examples are provided by models based on regime switching | 8 ], stochastic 
volatility [25] . jump processes [9], etc. Alternatively, it is possible to redistribute the risk by a temporally-fractional 
derivative, which, similar to spatial asymmetry, redistributes the risk to either the short-period options or long¬ 
term options. Such model enables one to treat the short-term differently; instant risk coming from contemporary 
fluctuations, jump corrections, etc., is treated differently from the long-term risk, which is determined by the slow 
volatility of the system. This also influences the optimal hedging policy. 

While the long-term average behavior of markets tends to be similar in systems driven by ordinary, first time- 
derivative diffusion (i.e., diffusion speed parameter equal to one), an adjustment of diffusion speed parameter 7 to 
values 7 7 ^ 1 can often describe the system better. Such fact coheres with the complex nature of financial markets and 
particularly option trading. Further investigations of the topic and its interrelation to other models will be the subject 
of future research. Connection with regime switching models or interpretation of 7 parameter as a “regime-switching” 
time-dependent parameter and the conjunction with other models are questions of high importance and can possibly 
reveal additional potential of models based on double-fractional diffusion. 
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Appendix A: Comparison of Smearing Kernels for Riesz-Feller and Caputo Derivatives 


In this appendix, we compare Green function of Double-fractional diffusion in case when time derivative is equal to 
Riesz-Feller derivative and Caputo derivative and 7 < 1. The Green function is given by 


-Ll.. ■ (— 


,00 ! 

9a n {^T) = j ' dl 

where is different according to used derivative, so 


^ Qol (^5 0 5 


1 \ J for Riesz derivative, 

M y for Caputo derivative. 


Indeed, it is interesting to see, what happens with the smearing kernel 


(Al) 


(A 2 ) 


(A3) 


for small and large values. Firstly, when l —)> 0, then the argument of 7 -stable distribution distribution goes to infinity. 
According to Ref. m , the asymptotic expansion gives us 


which gives us 


1 / t \ r( 7 + 1) sin(7T7) l 

l x h 7,1 \V-h) ~ cos (4^) r^ +1 


for l 0 , 




l r( 7 )T (7 + 1 ) sin( 7 T 7 ) 
^ cos ( 7 ) 


for l — 0 . 


For Caputo case we get non-zero value of smearing kernel for l = 0. Particularly, it is equal to 


(A4) 


(A5) 


9? (t, 0) 


f J_\ r(7) sin(7T7) 
W cos ( 7 ) 


(A 6 ) 


On the other hand, when l —)> 00 , then it is necessary to use the Taylor expansion of L 7 j i(x), and again from [43] we 
obtain 


L 7) i (x) ~ A 1 x 1 2 7 exp (— B 1 x A7 ) for x 0 + , (A7) 

where A 7 = and A 1 resp. B 1 are 7 -dependent constants. The asymptotic behavior can be therefore described as 

#f F ( r , 0 ~ C i?F (r)A 7 / 2 ( 1 -^) exp B 1 D{r)l^^j for l +00 (A8) 

and for Caputo case as 

g? ( T ? 0 ^ C c (r)A^l 2 ( 1 -7) _1 exp B 1 D(r)l T ^^ for / —>> 0 . (A9) 

The r-dependent constants C RF (r ), resp. C c (r) can be determined from previous expressions. Therefore, in both 
cases we become exponential decay in l. The graphs of both smearing kernels are depicted in Fig. [5] 








